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ABSTRACT 

In self-consistent iV-body simulations of collisionless systems, gravitational interac- 
tions are modified on small scales to remove singularities and simplify the task of 
numerically integrating the equations of motion. This 'gravitational softening' is some- 
times presented as an ad-hoc departure from Newtonian gravity. However, softening 
can also be described as a smoothing operation applied to the mass distribution; the 
gravitational potential and the smoothed density obey Poisson's equation precisely. 
While 'softening' and 'smoothing' are mathematically equivalent descriptions, the lat- 
ter has some advantages. For example, the smoothing description suggests a way to 
set up iV-body initial conditions in almost perfect dynamical equilibrium. 

Key words: methods: numerical - galaxies: kinematics & dynamics 



1 INTRODUCTION 

The evolution of a collisionless self-gravitating system is described by two coupled equations: the Vlasov equation, 

§£ +T .§£_v..g... (D 

where / = f(r,v,t) is the one-particle distribution function and $(r, t) is the gravitational potential, and Poisson's equation, 



V 2 $ = A-KGp = 4ttG / dv / 



(2) 

iV-body simula tions use a Monte-C arlo method to solve these equations. The distribution function is represented by a collection 
of N particles ( Klimontovicbl 1967 ) : 



> ; /(r,v,i) = Vm,i 3 (r-r,(i))i 3 (v-v,(t)) ; 



(3) 



where rm, r,, and are the mass, position, and velocity of particle i. Over time, particles move along characteristics of (TTJ); 
at each instant, their positions provide the density needed for ([2]). 

In many collisionless iV-body simulations, the equations of motion actually integrated are 



dvi 
~~dt 



~~dt 



N 



i |2 +£ 2)3/2 



(4) 



where e is the softening length. These equations reduce to the standard Newtonian equations of motion if e = 0. The main 
reason for setting e / is to suppress the 1 /r singularity in the Newtonian potential; this greatly simplifies the task of 



numerically integrating these equations (e.g., lDehnenll2001f h By limiting the spatial resolution of the gravitational force, 



softening also helps control fluctuations caused by sampling the dist ribution function with finite N\ however ; this comes at a 
price, since the gravitational field is systematically biased for e / ( Merritt 19961 : Athanassoula et al. 199&1 200Clh . 

Softening is often described as a modi fication of Newto nian gravity, with the 1/r potential replaced by l/\/r 2 + e 2 . The 
latter is proportional to the p otential of a Plummerl (|191ll ) sphere with scale radius e. This does not imply that particles 



interact like Plummer spheres ( Dyer fc Ip| 19931 ): the acceleration of particle i is comput ed from the field at the po int only 



But it does imply that softening can also be described as a smoothing operation (e.g., Hernquist fc Barnes! 1990h . in which 



the pointillistic Monte-Carlo representation of the density field is convolved with the kernel 
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Figure 1. Effect of Plummcr smoothing on a p oc r -1 profile. 
Dashed line is the underlying density profile; solid curve is 
the result of smoothing with e = 1. The smoothed profile is 
always less than the underlying power-law. 




Figure 2. Effect of Plummer smoothing on a p oc r~ 2 profile. 
Dashed line is the underlying density profile; solid curve is the 
result of smoothing with e = 1. The smoothed profile slightly 
exceeds the underlying power-law at large r. 



(5) 



(6) 



(7) 



S '( r '' e ) 47r ( r 2 + £ 2)5/2 ■ 

In effect, the source term for Poisson's equation Q is replaced with the smoothed density 
p(r;e) = / dr'p(r')S(|r - r'|; e) = J dv' p{v - r')5(|r'|; e) . 

Formally, @ provides a Monte-Carlo solution to the Vlasov equation JTJ coupled with 

V 2 $ = 47rGp(r; e) = 4ttG / dr J dv/(r',v,t)S(|r - r'|; e) . 

Thus one may argue that a softened iV-body simulation actually uses standard Newtonian gravity, as long as it is clear that 
the mass distribution generating the gravitational field is derived from the particles via a smoothing process. 

Although Plummer softening is widely used in iV-body simulations, its effects are incompletely understood. If the un- 
derlying density field is featureless on scales of order e, softening has relative ly little eff e ct. Ho weve r, iV-body simulations are 
often used to model systems with power-law density profiles; for example, Hernauist ( 199fj| ) and Navarro. Frenk. fc White! 
( 19961 'lereafter NFW) models, which have p oc r 1 at small r, are widely used as initial conditions. One purpose of this paper 
is to examine how softening modifies such profiles. 

Assume that the underlying density profile is spherically symmetric and centered on the origin: p = p(|r|). The integrand 
in ((6| is unchanged by rotation about the axis containing the origin and the point r, so the integral can be simplified by 
adopting cylindrical coordinates (<p, R, z), where r is located on the z-axis at z = |r|. The integral over ip is trivial; for Plummer 
smoothing, the result is 



p(r;e) 



2 



dz 



dRR 



p(VR 2 + z 2 ) 



(R 2 + (z - 



e 2)5/2 



3^ 
2 



dz 



dRR 



p(^/R2 + (z-r) 2 ) 
(R 2 + z 2 + e 2 fl 2 ' 



(8) 



where the second equality holds because the outer integral is taken over the entire z axis. 



2 POWER-LAW PROFILES 

The first step is to examine the effect of Plummer smoothing on power-law density profiles, p n (f) = p a (a/r) n , where < n < 3. 
These profiles are not realistic, since the total mass diverges as r — > oo. However, results obtained for power- law profiles help 
interpret the effects of smoothing on more realistic models. 



2.1 The case p oc r 1 

Let the density be pi(r) = p a (a/r). The total mass enclosed within radius r is Mi(r) = 2irp a ar 2 . In this case, the smoothed 
density profile can be calculated analytically; the double integral is 

, Jo (R 2 + z 2 + e 2 ) 6 / 2 3e 2 V?T^ 

This yields a remarkably simple result for the smoothed density, plotted in Fig. [TJ 
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Figure 3. Density ratio Do = pn(0; e)/pn(f) plotted as a function of n. Limiting values are Do = 1 as n — > and Do = oo as n — > 3. 

(10) 



pi(r;e) = p a 



Ve 2 + r 2 

where r e = \/e 2 + r 2 . The smoothed mass within radius r, hereafter called the smoothed mass profil^, is 



Mi(r;e) 



dxAnx 2 pi(x; e) = 27tp a a (Vr e — e 2 sinh 1 (r/e)) 



(11) 



2.2 The case p oc r 

Let the density be p2{r) = p a (a/r) 2 . The total mass enclosed within radius r is M2(r) = Aivp^r. The integral over R can be 
evaluated, but the result is not particularly informative and the remaining integral must be done numerically. Fig. [2] presents 
the results. For log(r) > 0.2, the smoothed density exceeds the underlying power-law profile. This occurs because smoothing, 
in effect, spreads mass from r < e to larger radii, and with the underlying profile dropping away so steeply this redistributed 
mass makes a relatively large contribution to p2(r; e). Note that as r — > 0, the smoothed density p2(r; e) — > 2p2(e). 



2.3 Central density 

It appears impossible to calculate the smoothed density profile for arbitrary n without resorting to numerical methods, but 
the central density is another matter. Setting r = 0, the smoothed density is 

pn(0;e) = 3e Vo dx (J+". 2 )5/ 2 =^ r (f-f)r(|K(E)- (12) 

The central density ratio Do(n) = pn{0; e)/p n (e) is plotted as a function of n in Fig. [3] For n = 1 and 2, the ratio Do = 1 and 
2, respectively, in accord with the results above, while as n -> 3 the central density diverges. 

The smoothed central density for an arbitrary power-law is useful in devising an approximate expression for the smoothed 
density profile (Appendix A.l). In addition, the central density is related to the shortest dynamical time-scale present in an 
JV-body simulation, which may in turn be used to estimate a maximum permissible value for the time-step (§ 4.3.1). 



3 ASTROPHYSICAL MODELS 
3.1 Hernquist and NFW models 



As no ted above, both of these profiles have p oc r 1 as r — > 0. For this reason, they are treated in parallel. The iHernquist 
(199(3) model has density and mass profiles 



/ \ aM „, . , Mr 2 , . 

P H ( r > = o — 7 — I — vT ' Mn(r) = - — ■ — ^ , (13) 
2-7rr(a + r) s (r + a) 2 



where a is the scale radius and M is the total mass. The lNavarro. Frenk. fc White! (|l996l ) model has density and mass profiles 



3 

PnfwM = , a P ° , MnfwM = 47rpoa 3 (log ( - — - ) — ) , (14) 

where a is again the scale radius and po is a characteristic density. The double integrals required to evaluate the smoothed 
versions of these profiles appear intractable analyticalljH but can readily be calculated numerically. Figs. [4] and \E\ present 

1 This profile can't be obtained by applying kernel smoothing directly to M(r); only density profiles can be smoothed. 

2 Smoothed central densities for these and other profiles can be expressed in terms of special functions. 
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Figure 4. Effect of Plummcr smoothing on Hcrnquist profile. 
Top curve shows the density profile of a Hcrnquist model with 
scale radius a = 1 and mass M = 1. Lower curves show profiles 
smoothed with e = 1/256, 1/128, . . . , 1 (from top to bottom); 
heavy curve is e = 1/64, dashed curve is e = 1. Inset shows 
ratio pn(r; e)/pu(r). 



Figure 5. Effect of Plummcr smoothing on NFW profile. Top 
curve shows the density profile of a NFW model with scale 
radius a = 1 and density po = 1/(2t). Lower curves show 
profiles smoothed with e = 1/256, 1/128, . . . , 1 (from top to 
bottom); heavy curve is e = 1/64, dashed curve is e = 1. Inset 
shows ratio PNFw( r ; £ )/pnfw(0- 



results for a range of e values between a and a/256. For comparison, both models are scaled to have the same underlying 
density profile at r < a. 

The smoothed profiles shown in Figs. [4] and [5] are, for the most part, easily understood in terms of the results obtained for 
power-laws. For radii r < e, smoothing transforms central cusps into constant-density cores, just as in Fig. [1] If the softening 
length e is much less than the scale length a, the smoothed density p(r; e) within r <C a is almost independent of the underlying 
profile at radii r > a. Consequently, the smoothed central density p(0; e) ~ p(e), echoing the result obtained for the power-law 
n = 1. In addition, the actual curves in Figs. [4] and [5] are shifted versions of the curves in Fig. [T] this observation motivates 
simple approximations to pn{r; e) and Pnfw(>~; e) described in Appendix A.2. 

On the other hand, if e is comparable to a, the quantitative agreement between these profiles and the smoothed n = 1 
profile breaks down; the smoothed density at small r has a non-negligible contribution from the underlying profile beyond the 
scale radius a. As an example, for e = 1 the central density of the smoothed NFW profile is higher than the central density 
of the smoothed Hernquist profile, because the former receives a larger contribution from mass beyond the scale radius. 

A somewhat more subtle result, shown in the insets, is that heavily smoothed profiles exceed the underlying profiles at 
radii r>a. This is basically the same effect found with the n = 2 power-law profile (§ 2.2); with the underlying density 
dropping rapidly as a function of r, the mass spread outward from smaller radii more than makes up for the mass spread to 
still larger radii. This effect is more evident for the Hernquist profile than for the NFW profile because the former falls off 
more steeply for r>a. 



3.2 Jaffe model 

The ljaffel |l983h model has density and mass profiles 

47rr 2 (a + r) A (r + a) 

where a is the scale radius and M is the total mass. The double integrals required to evaluate the smoothed version of this 
profile appear intractable analytically but can readily be calculated numerically. Fig. [6] present results for a range of e values 
between a and a/256. 

Again, much of the behavior shown in this plot can be understood by reference to the results for the n = 2 power-law. In 
particular, for smoothing lengths e«», the central density is pj(0; e) ~ 2pj(e), and the curves in Fig. [6] are shifted versions 
of the one in Fig. [2] As the inset shows, for larger values of e the smoothed profiles quite noticeably exceed the underlying 
profile; the effect is stronger here than it is for a Hernquist model because the Jaffe model has more mass within r < a to 
redistribute. 
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Figure 6. Effect of Plummcr smoothing on Jaffc profile. Top 
curve shows the density profile of a Jaffc model with scale 
radius a = 1 and mass M = 1. Lower curves show profiles 
smoothed with e = 1/256, 1/128, . . . , 1 (from top to bottom); 
heavy curve is e = 1/64, dashed curve is e = 1. Inset shows 
ratio pj(r;e)/pj(r). 




Figure 7. Logarithmic slopes of profiles from Figs. [4] [5] 
and [6] Solid lines are underlying profiles; from bottom to 
top, they represent Jaffe, Hernquist, and NFW models, re- 
spectively. Dashed, dot-dashed, and dotted lines give results 
for e = a/256, a/64, and a/16, respectively. 



3.3 How much softening is too much? 

Figs. [4] E] and[6]have interesting implications for iV-body experiments. One might expect the smoothed profiles to resolve the 
inner power-laws of the underlying models as long as the softening length e is somewhat less than the scale radius a, but that 
is not what these figures show. Profiles smoothed with e>a/16 are essentially constant-density cores attached to power-law 
outer profiles; the density within the core depends on e, but no inner cusp per se can be seen. For e<a/64, on the other hand, 
the smoothed profiles do appear to trace the inner power-laws over some finite range of radii, before flattening out at smaller 
r. Only for e< a/256 can the inner cusps be followed for at least a decade in radius. 

Fig. helps explain this result. The underlying Jaffe, Hernquist, and NFW profiles all roll over gradually from their inner 
to outer power-law slopes between radii 0.1o<r<10a. Thus a resolution somewhat better than 0.1a is required to see the 
inner cusps of these models. In practice, this implies the softening parameter e must be several times smaller than 0.1a. 



4 TESTS AND APPLICATIONS 

Since the formalism developed above is exact, numerical tests of a relation like ([8]) for the smoothed density p(r; e) may seem 
superfluous. In practice, such tests can be illuminating - as benchmarks of iV-body technique. In what follows, the smoothing 
formalism will be applied to actual iV-body calculations, to check JV-body methodology and to demonstrate that the formalism 
has real applications. 

Putting this plan into operation requires some care. To begin with, an JV-body realization of a standard Hernquist or Jaffe 
profile spans a huge range of radii. Typically, the innermost particle has radius ri n ~ aN -1 ^ 2 or aN" 1 for a Hernquist or Jaffe 
profile, respectively, while for either profile, the outermost particle has radius r out ~ aN . A dynamic range of 7W /fin ~ iV 3/2 
or N 2 can be awkward to handle numerically; even gridless tree codes may not accommodate such enormous ranges gracefully. 
One simple option is to truncate the particle distribution at some fairly large radius, but it's preferable to smoothly taper the 
density profile: 

f (l + /i)p(r), r < b 

p(r) -> pt(r) = i (16) 
{ (l + M)p*(oA) 2 e _r/r * , r>b 

where the taper radius b^> a, the values of r, and p* are fixed by requiring that pt(r) and its first derivative are continuous 
at r = b, and the value of p, <C 1 is chosen to preserve the total mass. Let 

P= T -f (17) 
p ar r=b 
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Figure 8. Effect of Plummcr smoothing on tapered Hernquist (left) and Jaffe (right) models, represented by the solid curves. Dashed, 
dot-dashed, and dotted curves show results for e = 1/256, 1/64, and 1/16, respectively. 



be the logarithmic slope of the density profile at r 
b 



p(b)e 



-(2+13) 



and 



b, and M(r) be the underlying mass profile; then 

M(oo) 



(18) 



— (2 + 13) ' ' ~™ n M(b) + 47r6 2 r,p(6) 

Fig. [5] shows how Plummer smoothing modifies tapered Hernquist and Jaffe profiles. Both profiles have scale radius 
= 1, taper radius b = 100, and mass M = 1; these parameters will be used in all subsequent calculations. In each case, the 
underlying profile follows the standard curve out to the taper radius b, and then rapidly falls away from the outer power law. 
At radii r < b, the smoothed profiles match those shown in Figs. [4] and [6l apart from the factor of (l + /i) used to preserve total 
mass. At larger radii, the smoothed profiles initially track the underlying tapered profiles, but then transition to asymptotic 
p oc r -5 power law tails. This occurs because the Plummer smoothing kernel ([5]) falls off as r~ 5 at large r; in fact, these power 
laws match p = 3Me 2 /4-7rr 5 , which is the large-r approximation for a point mass M smoothed with a Plummer kernel. The 
amount of mass in these r~ 5 tails is negligible. 



4.1 Gravitational potentials 

In principle, it's straightforward to verify that the smoothed profiles above generate potentials matching those obtained from 
iV-body calculations. For a given density profile p(r), construct a realization with particles at positions r;; a iV-body force 
calculation with softening e yields the gravitational potential $i for each particle. Conversely, given the smoothed density 
profile p(r; e), compute the smoothed mass profile M(r; e), and use the result to obtain the smoothed potential $(r; e): 

d$ _ M(r;e) 

dr ~ r 2 ' (19) 
with boundary condition <E> — > as r — > oo. For each particle, the iV-body potential $i may be compared with the predicted 
value $(|rj|;e); apart from V^V fluctuations, the two should agree. 

A major complication is that y/N fluctuations in iV-body realizations imprint spatially coherent perturbations on the 
gravitational field; potentials measured at adjacent positions are not statistically independent. For example, the softened 
potential at the origin of an iV-body system is 

Grrii 



(rf + e 2 )V2 



(20) 



If the radii rj are independently chosen, this expression is a Monte Carlo integral, which will deviate from $(0; e) by an amount 
of order 0(iV~ 1 ' /2 ); moreover, everywhere within r < e the potential will deviate upward or downward by roughly as much as it 
does at r = 0. One way around this is to average over many iV-body realizations, but this is tedious and expensive. An easier 
solution is to sample the radial distribution uniformly. Let M(r) be the mass profile associated with the underlying density 
p(r). Assign all particles equal masses, and determine the radius n of particle i by solving M(rj) = (i — 0.5)M(oo)/iV for 
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Figure 9. Difference <5<3?fc between A^-body and smoothed potentials for tapered Hernquist (left) and Jaffe (right) models, normalized 
by central potential <E>(0; e). Vertical dashed lines show value of e = 1/64. Points are results for uniform realizations; jagged curves shows 
averages for groups of 32 points. Grey-scale images show typical results for random realizations. Central potentials are $h(0; e) = —0.9809 
and *j(0;e) = -3.9009. 



i = 1 to N. This eliminates radial fluctuations; the Monte-Carlo integral for $o is replaced with a panel integration uniformly 
spaced in M(r), and the central potential is obtained with relatively high accuracy. 

This trick does not suppress non-radial fluctuations, so the TV-body potential evaluated at any point r > still differs 
from the true <3>(r; e). But a non-radial fluctuation which creates an overdensity at some position r OVC r must borrow mass from 
elsewhere on the sphere r = |r over |; over-dense and under-dense regions compensate each other when averaged over the entire 
surface of the sphere. The resulting potential fluctuations likewise average to zero over the sphere, as one can show by using 
Gauss's theorem to evaluate the average gradient of the potential and integrating inward from r = oo. 

Finally, a subtle bias arises if the particles used to generate the potential are also used to probe the potential, since 
local overdensities are sampled more heavily. To avoid this, the potential can be measured at a set of points r& which are 
independent of the particle positions r^. Then 5&k = &k — &{rk',e) should display some scatter, but average to zero when 
integrated over test points rj, within a spherical shell. 

Fig. [9] shows results from direct-sum potential calculations for tapered Hernquist and Jaffe models, using units with 
G — 1. In each case, the underlying density profile was represented with N = 2 20 = 1048576 equal-mass particles, and the 
potential was measured at 4096 points uniformly distributed in log(r) between O.Ole and lOOe. The points show results for 
uniform radial sampling. While non-radial fluctuations create scatter in <J?j,, the distribution is fairly symmetric about the 
line <5<3? = 0. The jagged curves are obtained by averaging 5$^ over radial bins each containing 32 points. These averages fall 
near zero, demonstrating very good agreement between the JV-body results and the potentials calculated from the smoothed 
density profiles. 

For comparison, the grey-scale images in Fig. [9] display representative results for random realizations of each density 
profile. In these realizations, the radius of particle i is computed by solving M(ri) = XiM{po), where Xi is a random numbeip 
uniformly distributed between and 1. To examine the range of outcomes, 1000 random realizations of each model were 
generated and ranked by central potential $ ; since particle radii are chosen independently, the central limit theorem implies 
that 4>o has a normal distribution. Shown here are the 25 th percentile and 75 th percentile members of these ensembles; half 
of all random realizations lie between the two examples presented in each figure. Note that these examples deviate from the 
true potential by fractional amounts of ~ A/~ 1//2 . Obviously, it's impossible to detect discrepancies between $fe and &(rk\ e) 
of less than one part in 10 3 using random realizations with N ~ 10 6 . 

It's instructive, not to mention disconcerting, to try reproducing Fig. [9] using a tree code ( Barnes fcFIut 19861) in- 
stead of direct summation. Tree codes employ approximations which become less accurate for e > ( HeTnauistl Il987 : 
Wachlin fc Carpintero! bOQfjh ; these systematically bias computed potentials and accelerations (see Appendix B). For ex- 



3 A good random number generator is essential. The Unix generator, randomO, appears to be slightly non-u niform; replacing X j with 
1 — Xi yields systematically different <E>o values. The results shown here use the Tausworthe generator taus2 l lGalassi et al.ll2009h . 



© 0000 RAS, MNRAS 000, 000-000 



8 J. E. Barnes 




-1 -0.8 -0.6 -0.4 -0.2 -6 -4 -2 

E E 

Figure 10. Distribution functions for tapered Hernquist (left) and Jaffe (right) models. Solid curve: /(£); dashed, dot-dashed, and 
dotted curves: f(E;e) for e = 1/256, 1/64, and 1/16, respectively. 



ample, the code which will shortly be used for dynamical tests, run with an opening angle 9 = 0.8, yields central potentials 
which are too deep by a few parts in 10 3 , depending on the system being modeled. This systematic error cannot be 'swept 
under the carpet' when comparing computed and predicted potentials at the level of precision attempted here. 



4.2 Distribution functions 

Constructing equilibrium configurations is an important element of many JV-body experiments. Approximate equilibria may 
be generated by a variety of ad hoc methods, but the construction of a true equilibrium iV-body model amounts to drawing 
particle positions and velocities from an equilibrium distribution function / = /(r, v). However, a configuration based on a 
distribution function (DF) derived without allowing for softening will not be in equilibrium if it is simulated with softening. 

Assume the model to be constructed is spherical and isotropic. Broadly speaking, there are two options: (a) adopt a DF 
/ = f(E) which depends on the ener gy _E ; and solve Poisson's equation for the gravitational potential, or (b) adopt a mass 



model p = p(r), and use Eddington's (|191g) formula to solve for the DF. If softening is taken into account, option (a) becomes 



somewhat awkward, since the source term for Poisson's equat ion (0 is non- locaB On the other hand, option (b) is relatively 
straightforward (e.g., Kazantzidis. Maeorrian. fc Moore! 20041 ). 



Starting with a desired density profile p(r), the first step is to compute the smoothed density and mass profiles p(r; e) 
and M(r;e), respectively. Since M(r; e) > everywhere, equation (|19|) guarantees that the smoothed potential <E>(r;e) is a 
monotonically increasing function of r. It is therefore possible to express the underlying density profile p(r) as a function of 
<£>(r; e), and compute the DF: 

m *'^Tgf, mi - Er1 '*^- (21) 

Note that in dp/d<&, the quantity p = p(r) is the underlying density, while <E> = <I>(r; e) is the smoothed potential, related by 
Poisson's equation to the smoothed density p(r;e). In effect, the smoothed potential $(r; e) is taken as a given, and (|21l) is 
used to find what will hereafter be called the s moothed distribution function f(E; e); with this DF, the underlying profile p(r) 
is in equilibrium in the adopted potential (e.g., McMillan fc Dehnenll2007 : Barnes fc Hibbard 20091 ) . Conversely, setting e = 



yields the self- consistent distribution function f(E) which describes a self-gravitating model with the underlying profile p(r). 

Fig.[l0]presents DFs f or tapered Hernquist and Ja ffe models. In each case the solid line shows the self-consistent DF; these 
match the published DFs (|jaffell 19831 : lHernciuisrfl990l ) over almost the entire energy range, deviating only for E> — GM/b = 



—0.01 where tapering sets in. 

Smoothed DFs for Jaffe models appear very different from their self-consistent counterpart. The latter has a logarithmic, 
infinitely deep potential well, which effectively confines material with constant velocity dispersion in a p oc r~ 2 cusp. The 
characteristic phase-space density / ~ po" 3 oc r~ 2 diverges as r — > (ie, as E — > -co), but only because p does. With e > 
the potential well is harmonic at small r, and can't confine a p oc r~ 2 cusp unless the local velocity dispersion scales as a oc r; 
thus the phase-space density now diverges as / oc r~ 5 . Moreover, the domain of f(E; e) is limited to to Eo < E < 0, where 
Eo = <E>(0; e). Thus, instead of growing exponentially as a function of —E, the smoothed DF abruptly diverges at some finite 
energy. 

By comparison, smoothed DFs for Hernquist models look similar to the self-consistent DF. The latter has a potential 
well of finite depth, and the smoothed profiles generate wells which are only slightly shallower. As the left panel of Fig. 1101 
shows, all the DFs asymptote to oo as E — > Eo. However, the run of velocity dispersion with r is different; the self-consistent 
model has a oc r 1//2 , implying / oc per -3 oc r~ 5 ^ 2 . In contrast, the smoothed models have a oc r, implying / oc r~ 



-4 



iDebattista fc Sellwoodl teOOCh describe an iterative scheme using softened Af-body potentials which implements option (a). 
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One consequence is that the way in which / — > oo as E — > Eq is different in the smoothed and self-consistent Hernquist 
models. The self-consistent model has a linear potential as small r, and thus / oc r _5//2 oc (E — Eo)~ 5 ^ 2 . By comparison, the 
models based on smoothed potentials have harmonic cores, and as a result, / oc r~ 4 oc (E — Ea)~ 2 . (This difference is not 
apparent in Fig. 1101 but becomes obvious when log(_B — Eo) is plotted against log /.) In this respect, the use of a smoothed 
potential effects a non-trivial change on Hernquist models: / is a different power-law of (E — Eq). Coincidentally, the smoothed 
Jaffe models have / oc r~ 5 oc (E — Eo)~ 5 ^ 2 , just like the self-consistent Hernquist model. 

4.3 Dynamical tests 

iV-body simulations are useful to show that the distribution functions just constructed are actually in dynamical equilibrium 
with their smoothed potentials. For each model and e value, two ensembles of three random realizations were run. In one 
ensemble, the initial conditions were generated using the self-consistent DF f(E). The other ensemble used initial conditions 
generated from the smoothed DF f(E; e), which allows for the effects of softening. 

Each realization contained TV = 2 18 = 262144 equal-mass particles. Initial particle radii n were selected randomly 
by solving M(vj) = XiM(od) as described above. Initial particle speeds m were selected randomly by rejection sampling 
( von Neumann|[l95ll ) from the distributions g(v; rj) = v 2 f(\v 2 + $(rj)) or g{v\ n) = v 2 f(\v 2 + < I > (r i ; e); e), where the former 
assumes the self-consistent DF, and the latter a smoothed DF. Position and velocity vectors for particle i are obtained by 
multiplying n and Vi by independent unit vectors drawn from an isotropic distribution. In effect, this procedure treats the 
6-D distribution function /(r, v) as a probability density, and selects each particle's coordinates independently. 

Simulations were run using a hierarchical TV-body codjf]. An opening angle of 9 = 0.8, together with quadrupole moments, 
provided forces with median errors <5a/|a| < 0.0006. Particles within r ~ lOe have much larger force errors, although these 
seem to have limited effect in practice (Appendix B). Trajectories were integrated using a time-centered leap-frog, with the 
same time-step At = 1/1024 for all particles (see § 4.3.1). All simulations were run to t = 16, which is more than sufficient to 
test initial equilibrium. 

Fig. [TT] shows how the potential well depth $ m i n of each simulation evolves as a function of time. Here, well depth is 
estimated by computing the softened gravitational potential <F; of each particle i and taking the minimum (most negative) 
value. (This is more accurate than evaluating the potential at the origin since the center of the system may wander slightly 
during a dynamical simulation.) To better display the observed changes in $ m m, the horizontal axis is logarithmic in time. 

Most of the ensembles set up without allowing for softening (dotted curves in Fig. Hip are clearly not in equilibrium. In 
all three of the Jaffe models (bottom row), the potential wells become dramatically shallower on a time-scale comparable to 
the dynamical time at r = e. The reason for this is evident. The self-consistent Jaffe model has a central potential diverging 
like logr as r — > 0; this potential can confine particles with finite velocity dispersion at arbitrarily small radii. However, the 
relatively shallow potential well of a smoothed Jaffe model cannot confine these particles; they travel outward in a coherent 
surge and phase-mix at radii of a few e. Their outward surge and subsequent fallback accounts for the rapid rise and partial 
rebound of the central potential. Similar although less pronounced evolution occurs in the Hernquist models (top row) with 
e = 1/16 and possibly with e = 1/64 as well. Only the self-consistent Hernquist models run with e = 1/256 appear truly close 
to equilibrium. 

In contrast, all of the ensembles set up with smoothed DFs (solid curves in Fig. Hip are close to dynamical equilibrium. 
In equilibrium, gravitational potentials fluctuate as individual particles move along their orbits. If particles are uncorrelated, 
the amplitude of these fluctuations should be comparable to the amplitude seen in an ensemble of independent realizations. 
To check this, 1000 realizations of each model were generated; central gravitational potentials $ m i n were evaluated using the 
same tree algorithm and parameters used for the self-consistent simulations. The grey horizontal bands show a range of ±2<r 
around the average central potential for each model and choice of e. Some of the simulations set up using smoothed DFs 
wander slightly beyond the 2a range. However, none of them exhibit the dramatic evolution seen in the cases set up using 
self-consistent DFs. 

The central potential is relatively insensitive to changes in the mass distribution on scales r <C e. To examine small-scale 
changes directly, density profiles measured from the initial conditions were compared to profiles measured at t = 16 time 
units. These profiles were derived as follows. First, SPH-style interpolation with an adaptive kernel containing 32 particles 
was used to estimate the density around each particle. Next, the centroid position r cont of the 32 highest-density particles was 
determined. Finally, a set of nested spherical shells, centered on r cen t, were constructed; each shell was required to contain at 
least 16 particles and have an outer radius at least 1.05 times its inner radius. 

Fig. H2] summarizes results for Jaffe models, which display the most obvious changes. The density of each shell is plotted 
against the average distance from r ccnt of the particles it contains. In each panel, the top set of curves compare initial 
(t — 0) numerical results with the underlying tapered Jaffe model, always represented by a light grey line. Profiles from three 

5 See http://www.ifa.hawaii.edu/faculty/barnes/treecode/treeguide.html for a description. 
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Figure 11. Evolution of potential well depth for iV-body simulations of Hcrnquist (top row) and Jaffe (bottom row) models, each run 
with the e value labeled. Solid (dotted) curves show results for initial conditions generated from smoothed (self-consistent) DFs. Light 
grey bands show expected ±2cr variation in central potential for N = 262144 independent particles. 
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Figure 12. Density profiles of Jaffe models before and after dynamical evolution. Initial conditions were constructed using self-consistent 
DFs (left) and smoothed DFs (right). In each panel, the top profile shows the initial conditions; the smooth curve is the underlying 
tapered model, overplottcd by three slightly bumpy curves showing numerical results from three independent realizations. The three 
profiles below show numerical results at t = 16 for simulations run with e = 1/256, 1/64, and 1/16, respectively; each is displaced 
downward by one additional unit in log p for clarity. 
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independent iV-body realizations of each model are overplotted. While some scatter from realization to realization is seen, the 
measured densities track the underlying profile throughout the entire range plotted. The outermost point is at ~ fO 7 times 
the radius of the innermost one; there are not enough particles to obtain measurements at smaller or larger radii. 

Ranged below the top curves in Fig. [12] are numerical results at t = 16 for softening lengths e = 1/256, 1/64, and 
1/16, each shifted downward by one more unit in log p. Again, profiles from three independent simulations are overplotted to 
illustrate run-to-run variations. Simulations set up using the self-consistent DF (left panel) show significant density evolution; 
their initial power- law profiles are rapidly replaced by cores of roughly constant density inside r <e. In contrast, simulations 
set up using smoothed DFs (right panel) follow the initial profile down to r ~ 10~ 4 (although density evolution occurs on 
smaller scales) . This shows that a careful set-up procedure can maintain the initial density profile even on scales much smaller 
than the softening length. 

A similar plot for the Hernquist models confirms that most of these simulations start close to equilibrium. Hernquist 
models set up using smoothed DFs don't appear to evolve at all, although this statement should be qualified since the profiles 
of these models can't be measured reliably on scales much smaller than r ~ 10~ 2 . Models set up using the self-consistent DF 
and run with e = 1/16 undergo some density evolution; their profiles fall below the underlying Hernquist model at r ~ e, 
although they continue rising to the innermost point measured. Simulations run with e = 1/64 or less display no obvious 
changes down to scales of r ~ 1CP 2 . 

It appears that the Jaffe models set up using smoothed DFs are not completely free of long-term evolution. The right- 
hand panel of Fig. [12] shows that the peak density as measured using a fixed number of particles falls by roughly an order 
of magnitude by t = 16. Moreover, a close examination of Fig. [11] turns some cases with a gradual decrease in potential well 
depth; in the smoothed Jaffe models with e = 1/256, for example, $ m i n exhibits an upward trend of ~ 0.011 percent per unit 
time. This evolution may not be due to any flaw in the initial conditions; the central cusps of such models, which are confined 
by very shallow harmonic potentials, are fragile and easily disrupted. There may be more than one mechanism at work here; 
the rate of potential evolution appears to be inversely proportional to particle number N, while the rate of density evolution 
is independent of N. A full examination of this matter is beyond the scope of this paper. 



4-3.1 Choice of time-step 

Selecting the time-step At for an iV-body simulation is a non-trivial problem. While the choice can usually be justified post-hoc 
by testing for convergence in a series of experiments with different time-steps, it's clearly convenient to be able to estimate 
an appropriate At a priori. A general rule governing such estimates is that the time-step should be smaller than the shortest 
dynamical time-scale present in the simulation. 

The central density p c = p(0; e) of a smoothed density profile defines one such time-scale. Within the nearly constant- 
density core of a smoothed profile, the local orbital period is t c = -^/37v/Gp c ; this is the shortest orbital period anywhere in 
the system. Numerical tests show the leap-frog integrator is well-behaved if At<0.05t c (conversely, it becomes unstable if 
At > 0.15t c )- Among the models simulated here, the Jaffe model with e = 1/256 has the highest smoothed central density; for 
this model, p c = pj(0;e) ~ D (2)pj(e) = 10440. Given this density, t c ~ 0.0300 and At < 0.0015. 

The time required for a fast-moving particle to cross the core region defines another time-scale. If $ c = $(0; e) is depth of 
the central potential well, the maximum speed of a bound particle is V— 2"3? c , and the core crossing time is t x = e/V— 2$T. The 
smoothed Jaffe model with e = 1/256 has the deepest potential well. For this model, tests of the leap-frog with fast-moving 
particles on radial orbits show that At < t x ~ 0.0012 yields good results, but time-steps a few times longer result in poor 
energy conservation as particles traverse the core region. (The relationship between tx and the maximum acceptable time-step 
At may be somewhat model-dependent.) 

Thus, for the Jaffe model with e = 1/256, both the local criterion based on t c and the global criterion based on t x 
yield similar constraint^] on At. It's convenient to round At down to the next power of two, implying At = 1/1024. This 
corresponds to At ~ 0.03t c — 0.81t x , which is somewhat conservative but helps insure that non-equilibrium changes will be 
accurately followed. 

To see if this time-step is reasonable, realizations of this model were simulated with various values of At between 1/128 
to 1/2048. At the lower end of this range, the effects of an over-large time-step manifest quickly; global energy conservation is 
violated, and the measured central potential $ m i n drifts upward over time (even though the initial conditions, generated from 
f(E;e), are near equilibrium). With a time-step At = 1/128, for example; the potential well becomes ~ 3 percent shallower 
during the first two time-steps, and by t — 4 its depth has decreased by 18 percent. These simulations also violate global 
energy conservation, becoming ~ 4.5 percent less bound by t = 4. Integration errors are reduced - but not entirely eliminated 
- with a time-step At = 1/256; by t — 4, the potential well becomes ~ 1.3 percent shallower, while total energy changes by 



6 Assuming e<a, one can show that t c /t x — ll^/ln(a/e) for any smoothed Jaffe model; both criteria yield similar constraints almost 
independent of e. For smoothed Hernquist models, on the other hand, t c /t x ~ 11a/ a/e; the constraint based on t x is generally stricter. 
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~ 0.4 percent. With a time-step of At = 1/512 or less, global energy conservation is essentially perfect, and variations in $ m i n 
appear to be driven largely by particle discreteness as opposed to time-step effects. 

Plots analogous to Fig. [T2]show that the simulations with time-steps as large as At = 1/256 ~ 0.12t c reproduce the inner 
cusps of Jaffe models just as well as those with At = 1/1024. With this time-step, individual particles may not be followed 
accurately, but their aggregate distribution is not obviously incorrect. On the other hand, a time-step At = 1/128 yields 
density profiles which fall below the initial curves for r <0.01. 



5 DISCUSSION 



Softening and smoothing are mathematically equivalent. While the particular form of softening adopted in ((4} corresponds to 
smo othing with a Plummer kernel |(5]), other softening prescriptions can also be described in terms of smoothing operations 
(e.g. lDehnenlboOll ). There are two conceptual advantages to thinking about softening as a smoothing process transforming the 



underlying density field p(r) to the smoothed density field p(r; e). First, since the gravitational potential <£>(r; e) is related to 
p(r; e) by Poisson's equation, the powerful mathematical machinery of classical potential theory becomes available to analyze 
potentials in iV-body simulations. Second, focusing attention on smoothing makes the source term for the gravitational field 
explicit. From this perspective, smoothing is not a property of the particles, but a separate step introduced to ameliorate 1/r 
singularities in the potential. Particle themselves are points rather than extended objects; this insures that their trajectories 
are characteristics of {TJ. 

Plummer smoothing converts p n (r) oc r~™ power-laws to cores. At radii r>e the density profile is essentially unchanged, 
while at r < e the density approaches a constant value equal to the density of the underlying model at r — e times a factor 
which depends only on n. For the case n = 1, this factor is unity and pi(r; e) = pi(\/r 2 + e 2 ) everywhere. 

The effects of Plummer smoothing on astrophysically-motivated models with power-law cusps, such as the Jaffe, Hernquist, 
and NFW profiles, follow for the most part from the results for pure power-laws. In particular, for e<a/64, where a is the 
profile's scale radius, the power-law results are essentially 'grafted' onto the underlying profile. On the other hand, for e > a/16, 
the inner power-law is erased by smoothing. 

Smoothing provides a way to predict the potentials obtained in iV-body calculations to an accuracy limited only by y/N 
fluctuations. These predictions offer new and powerful tests of iV-body methodology, exposing subtle systematic effects which 
may be difficult to diagnose by other means. 

Given an underlying density profile p(r), it's straightforward to construct an isotropic distribution function f(E;e) such 
that p(r) is in equilibrium with the potential generated by its smoothed counterpart p(r; e). Such distribution functions can be 
used to generate high-quality equilibri um initial conditions for iV-bo dy simulations; they should be particularly effective when 
realized with 'quiet start' procedures ( Debattista fc Sellwood 2000T I. Systems with shallow central cusps, such as Hernquist 



and NFW models, may be set up fairly close to equilibrium without taking softening into account as long as e is not too large. 
However, it appears impossible to set up a good iV-body realization of a Jaffe model without allowing for softening. 

It's true t hat realizations so constructed don't reproduce the actual dynamics of the underlying models at small radii 
( DehnenllioOl . botnote 8); to obtain an equilibrium, the velocity dispersion is reduced on scales r < e. But realizations set up 



without softening preserve neither the dispersion nor the density at small radii, and the initial relaxation of such a system 
can't be calculated a priori but must be simulated numerically. On the whole, it seems better to get the central density profile 
right on scales r < e, and know how the central velocity dispersion profile has been modified. Even if the dynamics are not 
believable within r < e, the ability to localize mass on such scales may be advantageous in modeling dynamics on larger scales. 

Mathematica code to tabulate smoothed models is available at 
http : //www. if a.hawaii . edu/f aculty/barnes/research/smoothing/. 
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APPENDIX A: APPROXIMATIONS 
A.l Power-law profiles 

The power-law density and cumulative mass profiles are 

Pn(r) = Pa (-)" , M n (r) = -i^p aa V 3 -> . (22) 
\r J 6 — n 

Plummer smoothing converts power laws with n < 3 to finite-density cores. At r < e the smoothed density is nearly 
constant and close to the smoothed central density p n {0; e) = Do(n)p n (e). Within this constant-density region, the smoothed 
mass profile is approximately 

Mn(r;e) = ^-r 3 p n (0;e) = ±^r 3 D (n)p a (J)" . (23) 

At r > e, on the other hand, smoothing has little effect on the mass profile, so M n (r; e) ~ M n (r). Interpolating between these 
functions yields an approximate expression for the smoothed mass profile: 

M n {r-e) = (M n (r;e)- K/ " + M n (r)- K/n p /K , (24) 

where the shape parameter k determines how abruptly the transition from one function to the other takes place. This expression 
can be rearranged to give 

^ ^ = ( ( (3-nUfro ) K ^ & K + 1 ) Mn(r) • (25) 
The smoothed density profile is obtained by differentiating the mass profile: 

Mr;e) = ^^M n (r;e). (26) 

Figs. [T31 and 1141 present tests of these approximations for p oc r^ 1 and r~ 2 power-laws, respectively. As in Figs. [1] and [2] 
the smoothed density profile was computed with e = 1; for other values of e, the entire pattern simply shifts left or right 
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Figure 13. Relative error in smoothed density (solid) and 
mass (dashed) for apar -1 profile, computed for e = 1 using 
H25H and (I26D . Dark curves show results for re = 1.739; light 
grey solid curves show errors in density only for re = 1.769 
(above) and ft = 1.709 (below). 
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Figure 14. Relative error in smoothed density (solid) and 
mass (dashed) for a p <x r~ 2 profile, computed for e = 1 using 
( 12 r >[ ) and H26I I. Dark curves show results for k = 1.820; light 
grey solid curves show errors in density only for n = 1.850 
(above) and n = 1.790 (below). 



without changing shape or amplitude. Dashed curves show relative errors in smoothed mass, Am = M (r; e) /M (r; e) — 1, 
while solid curves are relative errors in smoothed density A p = p(r; e)/p(r; e) — 1. The re value used for each dark curve is the 
value which minimizes y] . A p (ri) 2 evaluated at points n distributed uniformly in logr between logr = —1.5 and 1.5. In light 
grey, plots of A p (r) for two other re illustrate the sensitivity to this parameter. Comparing these plots, it appears that the 
approximation works better for the p oc r _1 power-law than it does for r -2 , but even in the latter case the maximum error is 
only ~ 2%. 

Because 1)251) modifies the underlying m ass profile with a multipl icative factor, it can also be used to approximate effects 
of softening on non-power-law profiles (e.g. Barnes fc Hibbard 20091 ) ; for this purpose, both e and re can be treated as free 
parameters and adjusted to provide a good fit. The resulting errors in density, which amount to a few percent near the 
softening scale, are undesirable but don't appear to seriously compromise A^-body simulations with N ~ 10 J . 



A. 2 Hernquist and NFW profiles 

For e <C a, smoothing primarily modifies the r _1 part of these density profiles. This, together with the exact solution for the 
case p oc r _1 given in § 2.1, suggests simple approximations for smoothed Hernquist and NFW models: 

pn(r; e) = pn(Vr 2 + e 2 ) = pn{r t ) , /5nfw(V; e) = /9nfw(\A" 2 + e 2 ) = pNFw(V e ) . (27) 

Fig. 1151 plots the relative error in density, A p = p(r; e) / p(r; e) — 1 for both models, adopting e = a/16. For other values of e, 
these errors scale roughly as e 1 ' 6 . 

The general behavior of these approximations is readily understood. Overall, pNFw( r ;e) is more accurate than pH(r;e) 
since the NFW profile is closer to p a r" 1 at all radii. Both curves are approximately flat for r <C e, then reach minima 
for r between e and the profile scale radius a. These minima arise because the smoothed density approaches or even slightly 
exceeds the underlying density (see Figs. [4] and [5}, while (|27|) always yields values below the underlying density. 

It's sometimes useful to have the cumulative mass for a smoothed profile. The approximate profiles in (|27[) can be 
integrated analytically, although the resulting expressions are a bit awkward: 

Af H (r;e) = / dx Atyx 1 pn{x; e) (28) 



arctan , — arctan 



rx / E 2_ a 2( a 3 T . £ _ a 2( e 2 +2r 2j + aT . e( . r 2_ e 2j + 1! 2 T .2^ 



and 



MNFw(r;e) = / da;47ra; 2 pNFw(2;;e) (29) 



+ ^ (arctan (_^=) - arctan (-^=) ) + log 



Note that because the approximate profiles (|27[) systematically underestimate the true smoothed densities, these expressions 
will likewise systematically underestimate the total smoothed mass. 
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Figure 15. Relative error in density, A p = p/p — 1, plotted vs. radius r, for the approximations given in I I27I I. Solid and dashed curves 
show results for Hcrnquist and NFW profiles, respectively, computed for e = a/16. 



APPENDIX B: FORCE CALCULATION ERRORS 



Tree codes reduce the computational cost of gravitational force calculation by making explicit approximations (jBarnes fc Hut 
19861 ). The long-range potential due to a localized mass distribution M with total mass M and center of mass position ro is 



approximated as 
*/ ^ GM 



+ higher order terms , 



(30) 



where the higher order terms include quadrupole and possibly higher-order moments (dipole terms vanish because ro coincides 
with the center of mass). To implement softening, this approximation is typically replaced with 
GM 



higher order terms . 



(31) 



y/\r - r () | 2 + e* 

This works at large dis tances, but becomes inaccura te if |r— ro| ~ e ( Hernquistl 1987 ) . Moreover, because the error is introduced 
at the monopole level ( Wachlin fc Carpintero 20061 ). higher-order corrections don't repair the damage. 

To appreciate the problem, consider a sphere 5 centered on ro with radius R large enough to enclose M. For e — 0, the 
inward acceleration averaged over the surface of 5 is easily computed using Gauss's theorem: 

If,. GM 



(32) 



In other words, the monopole term is sufficient to calculate the inward acceleration averaged over the surface of S exactly. 
Suppose we want to compute a r for e > 0. Again using Gauss's theorem, we have 

GM s {e) 



R? 



where Mj(e) = / drp(r;e) 



is the smoothed mass within the sphere. As before, this is an exact equality. The tree code approximation 
the enclosed mass is 
M 



M° s (e) = 



(1 + eVi? 2 ) 3 / 2 



(33) 

31) implies that 
(34) 



This is correct if M is simply a point mass located at ro. But if A4 has finite extent, then the enclosed mass Ms(e) is always 
less than M%(e). As a result. (I31C will systematically overestimate the inward acceleration and depth of the potential well 



due to M. Wachlin fc Carpintero ( 20061 ) demonstrate a similar result by computing the softened potential of a homogeneous 



sphere; they find — GM/yJ\r — ro| 2 + e 2 is only the first term in a series. 

The inequality Ms(e) < Mg(e) is easily verified for Plummer softening. An analogous inequality is lik ely to hold for other 
smoothing kernels S(r; e) which monotonically decrease with r. Smoothing kernels with compact support ( DehnenlbOQl ) may 
be better behaved in this regard. 

Under what conditions are these errors significant? For 'reasonable' values of e, most dynamically relevant interactions 
are on ranges Ar S> e where softening has little effect; these interactions are not compromised since (|31[) is nearly correct at 
long range. Only if a significant fraction of a system's mass lies within a region of size e can these err ors become important. 
This situation was not investigated in early tree code tests (e.g., Hernquist 1987 ; Barnes fc Hut 19891 ). which generally used 
mass models with cores instead of central cusps, and even heavily softened Hernquist models don't have much mass within 
one softening radius. On the other hand, Jaffe models pack more mass into small radii; a Jaffe model with e = a/16 has almost 
6 percent of its mass within r = e. Jaffe models should be good test configurations for examining treecode softening errors. 

Tests were run using tapered Jaffe and Hernquist models, realized using the same parameters (a = 1, b = 100, M = 1, 
and N = 262144) used in the dynamical experiments (§ 4.3). In each model, the gravitational field was sampled at 4096 points 
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Figure 16. Tree code acceleration error Sa/a plotted against 
radius. These results were obtained for a Jaffe model with 
9 = 0.8. Grey dots show errors for individual test points with 
e = 1/64; the jagged curve threading the dots is constructed 
by averaging points in groups of 16. Similar curves above and 
below show results for e = 1/16 and e = 1/256, respectively. 
The large marker on each curve shows the average value of 
Sa/a at radius r = e. The light grey curve shows results for 
e = 0. 




Figure 17. Jaffe model results for evolution of difference in 
potential well depth $ m i n for simulations run with 8 = 0.8 and 
0.4. Light grey, dark grey, and black show potential differences 
for e = 1/256, 1/64, and 1/16, respectively; three independent 
realizations arc plotted in each case. 



drawn from the same distribution as the mass. At each test point, results from the tree code with opening angle 9 = 0.8, 
including quadrupole terms, were compared with the results of an direct-sum code. As expected, the tests with Hernquist 
models showed relatively little trend of force calculation error with e, although the errors are somewhat larger for e = 1/16 
than for smaller values. In contrast, the tests with Jaffe models reveal a clear relationship between softening length and force 
calculation accuracy. 

Fig. 1161 shows the relative acceleration error Sa/a = |a t — ad|/|ad| for Jaffe models with various values of e. Here a t and 
ad are accelerations computed using a tree code and direct summation, respectively. The grey dots represent measurements 
of Sa/a for individual test points, computed using e = 1/64. The pattern of errors suggests two regimes. At radii r>0.25 
(log r > — 0.4), the points fall in a 'sawtooth' pattern which reflects the hierarchical cell structure used in the force calculation. 
At smaller radii, on the other hand, the relative error grows more or less monotonically as r — > 0. It appears that errors in 
the large-r regime are due to neglect of moments beyond quadrupole order in computing the potentials of individual cells; 
conversely, the errors in the small-r regime are due to the tree code's inaccurate treatment of softening. The direction of the 
error vectors a t — ad supports this interpretation; in the large-r regime they are isotropically distributed, while in the small-r 
regime they point toward the center of the system. 

The jagged line threading through the dots in Fig. 1161 constructed by averaging test points in groups of 16, shows the 
overall relationship between acceleration error and radius for e = 1/64. Similar curves are also plotted for e = 1/16 (above) 
and e = 1/256 (below). At large radii, all three curves coincide precisely, implying that force errors are independent of e. 
Going to smaller r, the curve for e = 1/16 is the first to diverge, rising above the other two, next the curve for e = 1/64 begins 
to rise, tracking the mean distribution of the plotted dots, and finally the curve for e = 1/256 parallels the other two. Each 
curve begins rising monotonically at a radius r ~ 20e; this is evidently where softening errors begin to dominate other errors 
in the force calculation. At the softening radius r = e, all three curves show mean acceleration errors Sa/a ~ 0.013. 

Are errors of this magnitude dynamically important? In particular, could they explain some of the potential evolution 
seen in the runs set up with softening (solid curves) in Fig. HIP One possible test is to re-run the simulations using smaller 
9 values; decreasing 9 from 0.8 to 0.6 or 0.4 reduces the tree code acceleration error associated with softening by factors of 
~ 2 or ~ 3, respectively (albeit at significant computational costs). The new runs were started using exactly the same initial 
conditions as their 9 — 0.8 counterparts. This initially allows central potentials $ m in(t) to be compared between runs to high 
accuracy, temporarily circumventing the effects of y~N fluctuations. 

Fig. 1171 compares Jaffe model results for 9 — 0.8 and 0.4. Initially, the potential difference $ m i n (t; = 0.8) — 4> m in(t; = 0.4) 
arises because the treecode systematically overestimates the potential well depth by an amount which is greater for larger 
values of 9. As the simulations run, the excess radial acceleration causes systems run with 9 = 0.8 to contract relative to 
those run with 9 — 0.4, causing a further increase in potential well depth. This contraction takes place on a dynamical 
time-scale at r ~ e, occurring first for the e = 1/256 simulations (light grey curves). At later times, as trajectories in 
otherwise-identical simulations with different 9 values diverge, short-term fluctuations in $ m m de-correlate and the dispersion 
in potential differences increases markedly. 

These results show that force-calculation errors can have measurable, albeit modest, effects on dynamical evolution. 
Extrapolating results for a range of 9 values to 6 = 0, it appears that for 9 = 0.8 the dynamical response of these Jaffe models 
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due to excess radial acceleration deepens central potential wells by about 0.2 percent. To test this extrapolation, it would 
be instructive to repeat these experiments with a direct-summation code. However, compared to the overall range of $ m i n 
variations seen in Fig. 1111 the perturbations due to force calculation errors seem relatively insignificant. In addition, density 
profiles measured from runs with 6 — 0.6 or 0.4 appear similar to those shown in Fig. 1121 again implying that treecode force 
calculation errors have little effect on the key results of this study. 

The best way to correct (|31[) for short-range interactions is not obvious. A simple, ad-hoc option is to reduce the effective 
opening angle on small scales. For example, accepting cells which satisfy 

d> -^-+5 where 6 cS = 9 p— c , (35) 

where d is the distance to the cell's center of mass , i is the cell' s size, measured along any edge, and S is the distance between 
the cell's center of mass and its geometric center ( Barnes 19981 ). reduces softening-relating treecode errors by a factor of ~ 2 
at a modest cost in computing time. Further experiments with similar expressions may produce better compromises between 
speed and accuracy. 
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